/*

   **********************************************************************
   *                                                                    *
   * This file contains the code for an image object for a content-     *
   * based image retrieval application that "summarizes" an image and   *
   * sends the resulting image "signature" to a server for              *
   * classification. The server returns the classification result in    *
   * the form of a Web page containing example images of the input      *
   * image's category.                                                  *
   *                                                                    *
   * The image class defined here is rather specific to this particular *
   * application and is not really designed for broad reuse (just yet). *
   *                                                                    *
   **********************************************************************
*/

/*
   **********************************************************************
   *                                                                    *
   * Gist   -- content-based image retrieval application for the gPhone *
   *                                                                    *
   * Copyright (C) 2011 Yong Zeng                                       *
   *                                                                    *
   * This file is part of Gist.                                       *
   *                                                                    *
   * Gist is free software; you can redistribute it and/or modify it  *
   * under the terms of the GNU General Public License as published by  *
   * the Free Software Foundation; either version 2 of the License, or  *
   * (at your option) any later version.                                *
   *                                                                    *
   * Gist is distributed in the hope that it will be useful, but      *
   * WITHOUT ANY WARRANTY; without even the implied warranty of         *
   * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU   *
   * General Public License for more details.                           *
   *                                                                    *
   * You should have received a copy of the GNU General Public License  *
   * along with Goggle; if not, write to the Free Software Foundation,  *
   * Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.       *
   *                                                                    *
   **********************************************************************
*/

//----------------------- PACKAGE SPECIFICATION -------------------------

package vision.gistcomputing ;

//------------------------------ IMPORTS --------------------------------

// Android utilities
import android.util.Log ;

// Java I/O support
import java.io.PrintStream ;
import java.io.FileOutputStream ;
import java.io.IOException ;

//------------------------- CLASS DEFINITION ----------------------------

/**
   This is a utility class for encapsulating the SIFT computation
   required by the BBoF estimator. Rather than implement the entire SIFT
   algorithm described in David Lowe's 2004 paper, this class goes
   straight to the SIFT descriptor part of it. That is, we skip the whole
   Gaussian scale space, key-point localization, etc. The Lazebnik
   algorithm does not require any of that stuff.
*/
class BbofSift {

/// Gaussians play a large role in the SIFT algorithm. This setting
/// specifies the standard deviation to use for these Gaussians.
private static final float SIGMA = 1.6f ;

/// SIFT descriptors are computed in pixel patches centered at a number
/// of key-points. This parameter specifies the half-size of this patch.
///
/// NOTE: This constant is called RADIUS because SIFT descriptor
/// computation also involves weighted gradient magnitudes that drop off
/// normally (i.e., Gaussian) from the key-point center in proportion to
/// the radial distance between a pixel in the patch and the center.
private static final float RADIUS = Math.round(5 * SIGMA) ;

/// The SIFT computations are applied on the input image supplied by the
/// BBoF estimator (this class's main and, for now at least, only)
/// client.
private Image m_input_image ;

//-------------------------- INITIALIZATION -----------------------------

/// Prior to applying SIFT on 16x16 pixel patches arranged in a regular
/// grid superimposed on the input image, we need to first remove all the
/// salt-n-pepper noise from the input image. Although Lazebnik, et al.
/// don't do this, we do because it seems like a good idea. The current
/// implementation of the BBoF estimator on the server side also first
/// denoises input images.
public BbofSift(Image I)
{
   m_input_image = Image.denoise(I) ;
}

//------------------------- SIFT COMPUTATIONS ---------------------------

/// This method applies the necessary steps of the SIFT algorithm so as
/// to obtain SIFT descriptors of 16x16 pixel patches of the input image.
/// Only the first scale and octave are used.
public Grid apply()
{
   Image.Gradient gradient = Image.gradient(m_input_image) ;

   final int W = m_input_image.width() ;
   final int H = m_input_image.height() ;
   final int R = (int) RADIUS ;

   Grid G = new Grid(W/R - 1, H/R - 1) ;
   for (int y = R; y < H - R + 1; y += R)
   for (int x = R; x < W - R + 1; x += R)
      G.add(sift_descriptor(x, y, gradient)) ;

   return G ;
}

// This method computes the SIFT descriptor of a 16x16 pixel patch
// centered at the point (X,Y).
private SiftDescriptor sift_descriptor(int X, int Y, Image.Gradient G)
{
   // Just some precomputed constants that allow us to minimize floating
   // point operations inside this function's inner loop.
   final float _1_OVER_2SIGMA_SQR = -0.5f/(RADIUS * RADIUS) ;
   final float _2PI               = (float)(2 * Math.PI) ;
   final float _1_OVER_PI_BY_4    = (float)(4/Math.PI) ;
   final float _2_OVER_R          = 2/RADIUS ;

   final int W = m_input_image.width() ;
   final int H = m_input_image.height() ;
   final int R = (int) RADIUS ;

   FeatureVector feature_vector = new FeatureVector() ;

   // Iterate over a square pixel patch of size 2*R + 1 centered at (X,Y)
   // and extract the SIFT features in this area of the image by binning
   // the gradients therein and weighting the bins using a Gaussian that
   // drops off according to the radial distance between each point in
   // the patch and its center (X,Y).
   for (int y = -R; y < R; ++y)
   for (int x = -R; x < R; ++x)
   {
      final int img_x = X + x ;
      final int img_y = Y + y ;
      if (img_x < 0 || img_x >= W || img_y < 0 || img_y >= H)
         continue ; // image coordinates out of bounds; try next pixel...

      // Find the dominant orientation at (x,y). This is the SIFT feature
      // for this pixel. This orientation is assigned to one (or more) of
      // eight possible bins (i.e., 0-45, 45-90, etc.).
      //
      // NOTE: All the range checking with PI is just to bring the angle
      // into the [0, 360] range.
      float gradient_angle = G.direction.get(img_x, img_y) % _2PI ;
      if (gradient_angle < 0)
         gradient_angle += _2PI ;
      if (gradient_angle >= Math.PI)
         gradient_angle -= _2PI ;
      final float orientation = // dominant orientation split into 8 bins
         4 + _1_OVER_PI_BY_4 * gradient_angle ;

      // The feature at (x,y) will be normally (i.e., Gaussian) weighted
      // based on the distance between (x,y) and (X,Y).
      final float exp = (float) Math.exp((x*x + y*y) * _1_OVER_2SIGMA_SQR) ;
      final float weight = exp * G.magnitude.get(img_x, img_y) ;

      // Add the feature to the feature vector at fractional coordinates
      // corresponding to a 4x4 grid of bins.
      final float xf = 2 + _2_OVER_R * x ;
      final float yf = 2 + _2_OVER_R * y;
      feature_vector.add(xf, yf, orientation, weight) ;
   }
   return new SiftDescriptor(feature_vector.to_sift()) ;
}

//----------------------- SIFT FEATURE VECTORS --------------------------

/**
   This class acts as a temporary store for SIFT key-point features. The
   final feature vector in a SIFT key-point is an array of bytes.
   However, here, we use floats to allow normalizations and other such
   operations.

   As detailed in David Lowe's 2004 SIFT paper, a SIFT descriptor is
   essentially a histogram with 128 bins (4 bins for x * 4 bins for y * 8
   bins for orientation).
*/
private class FeatureVector {

/// The sizes of the different bins.
private static final int X_SIZE = 4 ;
private static final int Y_SIZE = 4 ;
private static final int Z_SIZE = 8 ;

// Some precomputed values to try and speed things up.
private static final float X_SIZE_BOUNDARY = X_SIZE - 0.5f ;
private static final float Y_SIZE_BOUNDARY = Y_SIZE - 0.5f ;
private static final float Z_SIZE_BOUNDARY = Z_SIZE - 0.5f ;

private static final int ZY_SIZE = Z_SIZE * Y_SIZE ;

/// As detailed in David Lowe's 2004 SIFT paper, the descriptor values
/// are normalized, thresholded and then renormalized. This parameter
/// specifies the threshold limit to use.
private static final float THRESHOLD = 0.2f ;

/// The SIFT features are stored in an array of floats whose size equals
/// the total number of bins.
private float[] m_features = new float[X_SIZE * Y_SIZE * Z_SIZE] ;

/// This method adds the specified value to the appropriate bins of the
/// SIFT descriptor.
public void add(float x, float y, float z, float value)
{
   int   xi0, xi1, yi0, yi1, zi0, zi1 ; // the bins
   float wx0, wy0, wz0, wx1, wy1, wz1 ; // the corresponding weights

   // If close to bounds then the values go fully into the end bins,
   // otherwise they split between two adjacent bins. Note: a value of
   // 2.0 should equally split between bins 1 and 2.
   if (x <= 0.5f)
   {
      xi0 = 0 ;
      xi1 = 0 ;
      wx0 = 0.5f ;
      wx1 = 0.5f ;
   }
   else if (x >= X_SIZE_BOUNDARY)
   {
      xi0 = X_SIZE - 1 ;
      xi1 = X_SIZE - 1 ;
      wx0 = 0.5f ;
      wx1 = 0.5f ;
   }
   else
   {
      final float xx = x - 0.5f ;
      xi0 = (int) xx ;
      xi1 =  1 + xi0 ;
      wx1 = xx - xi0 ;
      wx0 =  1 - wx1 ;
   }
   if (y <= 0.5f)
   {
      yi0 = 0 ;
      yi1 = 0 ;
      wy0 = 0.5f ;
      wy1 = 0.5f ;
   }
   else if (y >= Y_SIZE_BOUNDARY)
   {
      yi0 = Y_SIZE - 1 ;
      yi1 = Y_SIZE - 1 ;
      wy0 = 0.5f ;
      wy1 = 0.5f ;
   }
   else
   {
      final float yy = y - 0.5f ;
      yi0 = (int) yy ;
      yi1 =  1 + yi0 ;
      wy1 = yy - yi0 ;
      wy0 =  1 - wy1 ;
   }
   // The situation is slightly different for orientation (i.e., z)
   // because we wrap around.
   if (z <= 0.5f)
   {
      zi0 = 0 ;
      zi1 = Z_SIZE - 1 ;
      wz0 =   0.5f + z ;
      wz1 =   1.0f - wz0 ;
   }
   else if (z >= Z_SIZE_BOUNDARY)
   {
      zi0 = Z_SIZE - 1 ;
      zi1 = 0 ;
      wz0 = Z_SIZE + 0.5f - z ;
      wz1 = 1 - wz0 ;
   }
   else
   {
      final float zz = z - 0.5f ;
      zi0 = (int) zz ;
      zi1 =  1 + zi0 ;
      wz1 = zz - zi0 ;
      wz0 =  1 - wz1 ;
   }

   // Convention: we add 1 for each unit of orientation (our fastest
   // varying index), then zSize for each unit of y, finally zSize*ySize
   // for each unit of x.
   xi0 *= ZY_SIZE ;
   xi1 *= ZY_SIZE ;
   yi0 *= Z_SIZE ;
   yi1 *= Z_SIZE ;

   m_features[xi0 + yi0 + zi0] += value * wx0 * wy0 * wz0 ;
   m_features[xi1 + yi0 + zi0] += value * wx1 * wy0 * wz0 ;
   m_features[xi0 + yi1 + zi0] += value * wx0 * wy1 * wz0 ;
   m_features[xi1 + yi1 + zi0] += value * wx1 * wy1 * wz0 ;
   m_features[xi0 + yi0 + zi1] += value * wx0 * wy0 * wz1 ;
   m_features[xi1 + yi0 + zi1] += value * wx1 * wy0 * wz1 ;
   m_features[xi0 + yi1 + zi1] += value * wx0 * wy1 * wz1 ;
   m_features[xi1 + yi1 + zi1] += value * wx1 * wy1 * wz1 ;
}

/// This method converts the SIFT features computed by add() to the final
/// byte values that will make up the SIFT descriptor for each pixel
/// patch. The features are first normalized, thresholded and then
/// renormalized prior to the conversion.
public int[] to_sift()
{
   Utils.normalize2(m_features) ;
   Utils.threshold(m_features, THRESHOLD) ;
   Utils.normalize2(m_features) ;

   int[] sift = new int[m_features.length] ;
   for (int i = 0; i < sift.length; ++i)
      sift[i] = Utils.clamp_byte(m_features[i] * 512) ;
   return sift ;
}

} // end of inner class BbofSift.FeatureVector

//----------------------------- SIFT GRID -------------------------------

/// This class encapsulates the functionality of a grid of SIFT
/// descriptors.
public class Grid {

/// The grid is implemented as a two-dimensional array.
private SiftDescriptor[][] m_grid ;

/// These indices together act as an iterator for adding new descriptors
/// to the grid.
private int x, y ;

/// Create an empty grid of the specified dimensions.
public Grid(int W, int H)
{
   m_grid = new SiftDescriptor[H][W] ;
   x = y = 0 ;
}

/// Add a new SIFT descriptor to the next available slot in the grid.
public void add(SiftDescriptor S)
{
   if (space_available())
   {
      m_grid[y][x] = S ;

      ++x ;
      if (x >= m_grid[0].length) {
         x = 0 ;
         ++y ;
      }
   }
   else
      throw new ArrayIndexOutOfBoundsException("SIFT grid is full") ;
}

// Return true if the iterator indices indicate that the grid is not yet
// full.
private boolean space_available()
{
   return x >= 0 && x < m_grid[0].length
       && y >= 0 && y < m_grid.length ;
}

/// Retrieve the width of the SIFT grid.
public int width()
{
   return m_grid[0].length ;
}

/// Retrieve the height of the SIFT grid.
public int height()
{
   return m_grid.length ;
}

/// Retrieve the SIFT descriptor at position (x,y).
public SiftDescriptor get(int x, int y)
{
   return m_grid[y][x] ;
}

/// Debug support (dumps the first couple of SIFT descriptors in the
/// grid.
public void dump()
{
   Log.e("Goggle", "Grid.dump(): size = " + width() + "x" + height()) ;
   final int n = Math.min(width(), height()) ;
   for (int i = 0; i < n; ++i)
      get(i, i).dump() ;
}

/// Debug support: save the SIFT grid to a plain text file, which can be
/// pulled from the phone and compared to some ground truth.
public void save(String file_name)
{
   try
   {
      PrintStream out = new PrintStream(new FileOutputStream(file_name)) ;
      for (int y = 0; y < height(); ++y)
      for (int x = 0; x < width();  ++x)
         out.println(get(x, y).toString()) ;
      out.close() ;
   }
   catch (IOException e)
   {
      Log.e("Goggle", "unable to save float array to " + file_name, e) ;
   }
}

} // end of inner class BbofSift.Grid

//-----------------------------------------------------------------------

} 
